subroutine beta_alkene_init( moltype )
use prec
use cell
use general
implicit none
   !----------------------------------------------------------------------------
   integer, intent(in) :: moltype
   !----------------------------------------------------------------------------
   integer  :: i
   !----------------------------------------------------------------------------
   subnamep = subname
   subname  = 'beta_alkene_init'
   !
   if ( allocated( bondsMol  ) ) deallocate( bondsMol  )
   if ( allocated( anglsMol  ) ) deallocate( anglsMol  )
   if ( allocated( dihrsMol  ) ) deallocate( dihrsMol  )
   if ( allocated( attypMol  ) ) deallocate( attypMol  )
   if ( allocated( atposMol  ) ) deallocate( atposMol  )
   if ( allocated( molIDMol  ) ) deallocate( molIDMol  )
   if ( allocated( rotatedMol) ) deallocate( rotatedMol)
   if ( allocated( chargeMol ) ) deallocate( chargeMol )
   if ( allocated( EMass     ) ) deallocate( EMass     )
   if ( allocated( sigma     ) ) deallocate( sigma     )
   if ( allocated( epison    ) ) deallocate( epison    )
   if ( allocated( bondcoef  ) ) deallocate( bondcoef  )
   if ( allocated( anglcoef  ) ) deallocate( anglcoef  )
   if ( allocated( torscoef  ) ) deallocate( torscoef  )
   if ( allocated( DistMin   ) ) deallocate( DistMin   )
   if ( allocated( EName     ) ) deallocate( EName     )
   !
   ntypeMol   = 4  ! 1: CT, Ct; 2: CM, Cm; 3: HC, Hc; 4: HT
   nbndtypMol = 5  ! 1: CT-CT; 2: CT-HT; 3: CM-CT; 4: CM-CM; 5: CM-HC
   nangtypMol = 8  ! 1: CT-CT-CT; 2: CT-CT-HT; 3: HC-CT-HC;
                   ! 4: CM-CT-CT; 5: CM-CT-HC; 6: CM-CM-CT
                   ! 7: CM-CM-HC; 8: HC-CM-HC, CT-CM-HC
   ndihtypMol = 7  ! 1: CT-CT-CT-CT, CT-CT-CT-CM; 2: CT-CT-CT-HC, HC-CT-CT-HC
                   ! 3: HC-CT-CT-CM; 4: HC-CM-CT-CT, HC-CM-CT-HC
                   ! 5: CM-CM-CT-CT; 6: HC-CT-CM-CM; 7: HC-CM-CM-HC, HC-CM-CM-CT
   !
   allocate( EMass(ntypeMol), sigma(ntypeMol), epison(ntypeMol),           &
   &         bondcoef(2, nbndtypMol), anglcoef(2, nangtypMol),             &
   &         torscoef(4, ndihtypMol), DistMin(ntypeMol), EName(ntypeMol+1) )
   EMass(1:2) = 12.0107D0
   EMass(3:4) = 1.00794D0
   sigma      = (/ 3.50D0,  3.55D0,  2.42D0, 2.50D0 /)
   epison     = (/ 0.066D0, 0.076D0, 0.030D0, 0.030D0 /)
   DistMin    = (/ 1.7D0, 1.7D0, 1.1D0, 1.1D0 /)
   EName      = (/ 'C ', 'C ', 'H ', 'H ', 'Cu' /)
   !
   bondcoef = reshape( (/ 268.D0, 1.529D0,  & ! 1 CT-CT
   &                      340.D0, 1.09D0,   & ! 2 CT-HC
   &                      317.D0, 1.51D0,   & ! 3 CM-CT
   &                      549.D0, 1.34D0,   & ! 4 CM-CM
   &                      340.D0, 1.08D0/), & ! 5 CM-HC
   &                     (/2, nbndtypMol/) )
   anglcoef = reshape( (/ 58.35D0,112.7D0,  & ! 1 CT-CT-CT
   &                      37.5D0, 110.7D0,  & ! 2 CT-CT-HC
   &                      33.D0,  107.8D0,  & ! 3 HC-CT-HC
   &                      63.0D0, 111.1D0,  & ! 4 CM-CT-CT
   &                      35.D0,  109.5D0,  & ! 5 CM-CT-HC
   &                      70.D0,  124.0D0,  & ! 6 CM-CM-CT
   &                      35.D0,  120.0D0,  & ! 7 CM-CM-HC
   &                      35.D0,  117.0D0/),& ! 8 HC-CM-HC
   &                      (/2, nangtypMol/) )
   torscoef = reshape( (/ 1.300D0, -0.050D0,  0.200D0,  0.0D0, & ! 1 CT-CT-CT-CT CT-CT-CT-CM
   &                      0.000D0,  0.000D0,  0.300D0,  0.0D0, & ! 2 HC-CT-CT-CT HC-CT-CT-HC
   &                      0.000D0,  0.000D0,  0.366D0,  0.0D0, & ! 3 CM-CT-CT-HC
   &                      0.000D0,  0.000D0,  0.318D0,  0.0D0, & ! 4 HC-CM-CT-HC HC-CM-CT-CT
   &                      0.346D0,  0.405D0, -0.904D0,  0.0D0, & ! 5 CM-CM-CT-CT
   &                      0.0D0,    0.0D0,   -0.372D0,  0.0D0, & ! 6 HC-CT-CM-CM
   &                      0.0D0,   14.0D0,    0.0D0,    0.0D0/)& ! 7 CT-CM-CM-HC HC-CM-CM-CT
   &                      , (/4, ndihtypMol/) )
   !
   select case ( moltype )
   case ( 1 ) ! beta hexene
      molform    = 'H3C-CH=CH-CH2-CH2-CH3'
      molname    = 'Hexene_beta'
      den0       = 0.673D0
      natomMol   = 18
      nbondsMol  = 17
      nanglsMol  = 30
      ndihrsMol  = 34
      !
      allocate( bondsMol(3, nbondsMol), AnglsMol(4, nanglsMol),                &
      &         dihrsMol(5, ndihrsMol), attypMol(natomMol),                    &
      &         atposMol(3,natomMol), chargeMol(natomMol), molIDMol(natomMol) )                
      !
      attypMol = (/ 1, 4, 4, 4, 2, 3, 2, 3, 1, 4, 4, 1, 4, 4, 1, 4, 4, 4 /)
      atposMol(:, 1) = (/ -8.8443D0,   0.1349D0,  0.0000D0 /) ! C --> H3C
      atposMol(:, 2) = (/ -9.6813D0,   0.8331D0,  0.0000D0 /) ! H --> H3C
      atposMol(:, 3) = (/ -8.8952D0,  -0.4925D0,  0.8899D0 /) ! H --> H3C
      atposMol(:, 4) = (/ -8.8952D0,  -0.4925D0, -0.8899D0 /) ! H --> H3C
      atposMol(:, 5) = (/ -7.5197D0,   0.9204D0,  0.0000D0 /) ! C --> =CH
      atposMol(:, 6) = (/ -7.5324D0,   2.0103D0,  0.0000D0 /) ! H --> =CH
      atposMol(:, 7) = (/ -6.1772D0,   0.1660D0,  0.0000D0 /) ! C --> =CH
      atposMol(:, 8) = (/ -6.1645D0,  -0.9240D0,  0.0000D0 /) ! H --> =CH
      atposMol(:, 9) = (/ -4.8526D0,   0.9514D0,  0.0000D0 /) ! C --> CH2 4
      atposMol(:,10) = (/ -4.8599D0,   1.5805D0, -0.8901D0 /) ! H --> CH2 4
      atposMol(:,11) = (/ -4.8599D0,   1.5805D0,  0.8901D0 /) ! H --> CH2 4
      atposMol(:,12) = (/ -3.5101D0,   0.1970D0,  0.0000D0 /) ! C --> CH2 5
      atposMol(:,13) = (/ -3.5028D0,  -0.4321D0,  0.8901D0 /) ! H --> CH2 5
      atposMol(:,14) = (/ -3.5028D0,  -0.4321D0, -0.8901D0 /) ! H --> CH2 5
      atposMol(:,15) = (/ -2.1855D0,   0.9824D0,  0.0000D0 /) ! C --> CH3 6
      atposMol(:,16) = (/ -1.3485D0,   0.2843D0,  0.0000D0 /) ! H --> CH3 6
      atposMol(:,17) = (/ -2.1346D0,   1.6099D0,  0.8899D0 /) ! H --> CH3 6
      atposMol(:,18) = (/ -2.1346D0,   1.6099D0, -0.8899D0 /) ! H --> CH3 6
      !
      chargeMol = (/ -0.180D0, 0.060D0,  0.060D0,  0.060D0, -0.115D0,  &
      &              0.115D0, -0.115D0,  0.115D0, -0.120D0,  0.060D0,  &
      &              0.060D0, -0.120D0,  0.060D0,  0.060D0, -0.180D0,  &
      &              0.060D0,  0.060D0,  0.060D0                      /) 
      !                  Type,  atom 1, atom 2
      bondsMol(:, 1) = (/ 2,  2,  1 /)
      bondsMol(:, 2) = (/ 2,  3,  1 /)
      bondsMol(:, 3) = (/ 2,  4,  1 /)
      bondsMol(:, 4) = (/ 3,  5,  1 /)
      bondsMol(:, 5) = (/ 5,  6,  5 /)
      bondsMol(:, 6) = (/ 4,  7,  5 /)
      bondsMol(:, 7) = (/ 5,  8,  7 /)
      bondsMol(:, 8) = (/ 3,  9,  7 /)
      bondsMol(:, 9) = (/ 2, 10,  9 /)
      bondsMol(:,10) = (/ 2, 11,  9 /)
      bondsMol(:,11) = (/ 1, 12,  9 /)
      bondsMol(:,12) = (/ 2, 13, 12 /)
      bondsMol(:,13) = (/ 2, 14, 12 /)
      bondsMol(:,14) = (/ 1, 15, 12 /)
      bondsMol(:,15) = (/ 2, 16, 15 /)
      bondsMol(:,16) = (/ 2, 17, 15 /)
      bondsMol(:,17) = (/ 2, 18, 15 /)
      !                   Type,  atom1,  atom2, atom3
      AnglsMol(:, 1) = (/ 3,  2,  1,  3 /)
      AnglsMol(:, 2) = (/ 3,  2,  1,  4 /)
      AnglsMol(:, 3) = (/ 5,  2,  1,  5 /)
      AnglsMol(:, 4) = (/ 3,  3,  1,  4 /)
      AnglsMol(:, 5) = (/ 5,  3,  1,  5 /)
      AnglsMol(:, 6) = (/ 5,  4,  1,  5 /)
      AnglsMol(:, 7) = (/ 8,  1,  5,  6 /)
      AnglsMol(:, 8) = (/ 6,  1,  5,  7 /)
      AnglsMol(:, 9) = (/ 7,  6,  5,  7 /)
      AnglsMol(:,10) = (/ 7,  5,  7,  8 /)
      AnglsMol(:,11) = (/ 6,  5,  7,  9 /)
      AnglsMol(:,12) = (/ 8,  8,  7,  9 /)
      AnglsMol(:,13) = (/ 5,  7,  9, 10 /)
      AnglsMol(:,14) = (/ 5,  7,  9, 11 /)
      AnglsMol(:,15) = (/ 4,  7,  9, 12 /)
      AnglsMol(:,16) = (/ 3, 10,  9, 11 /)
      AnglsMol(:,17) = (/ 2, 10,  9, 12 /)
      AnglsMol(:,18) = (/ 2, 11,  9, 12 /)
      AnglsMol(:,19) = (/ 2,  9, 12, 13 /)
      AnglsMol(:,20) = (/ 2,  9, 12, 14 /)
      AnglsMol(:,21) = (/ 1,  9, 12, 15 /)
      AnglsMol(:,22) = (/ 3, 13, 12, 14 /)
      AnglsMol(:,23) = (/ 2, 13, 12, 15 /)
      AnglsMol(:,24) = (/ 2, 14, 12, 15 /)
      AnglsMol(:,25) = (/ 2, 12, 15, 16 /)
      AnglsMol(:,26) = (/ 2, 12, 15, 17 /)
      AnglsMol(:,27) = (/ 2, 12, 15, 18 /)
      AnglsMol(:,28) = (/ 3, 16, 15, 17 /)
      AnglsMol(:,29) = (/ 3, 16, 15, 18 /)
      AnglsMol(:,30) = (/ 3, 17, 15, 18 /)
      !                      Type, atom1, atom2, atom3, atom4
      DihrsMol(:, 1) = (/ 4,  6,  5,  1,  2 /)
      DihrsMol(:, 2) = (/ 4,  6,  5,  1,  3 /)
      DihrsMol(:, 3) = (/ 4,  6,  5,  1,  4 /)
      DihrsMol(:, 4) = (/ 6,  7,  5,  1,  2 /)
      DihrsMol(:, 5) = (/ 6,  7,  5,  1,  3 /)
      DihrsMol(:, 6) = (/ 6,  7,  5,  1,  4 /)
      DihrsMol(:, 7) = (/ 7,  8,  7,  5,  1 /)
      DihrsMol(:, 8) = (/ 7,  8,  7,  5,  6 /)
      DihrsMol(:, 9) = (/ 7,  9,  7,  5,  1 /)
      DihrsMol(:,10) = (/ 7,  9,  7,  5,  6 /)
      DihrsMol(:,11) = (/ 6, 10,  9,  7,  5 /)
      DihrsMol(:,12) = (/ 4, 10,  9,  7,  8 /)
      DihrsMol(:,13) = (/ 6, 11,  9,  7,  5 /)
      DihrsMol(:,14) = (/ 4, 11,  9,  7,  8 /)
      DihrsMol(:,15) = (/ 5, 12,  9,  7,  5 /)
      DihrsMol(:,16) = (/ 4, 12,  9,  7,  8 /)
      DihrsMol(:,17) = (/ 3, 13, 12,  9,  7 /)
      DihrsMol(:,18) = (/ 2, 13, 12,  9, 10 /)
      DihrsMol(:,19) = (/ 2, 13, 12,  9, 11 /)
      DihrsMol(:,20) = (/ 3, 14, 12,  9,  7 /)
      DihrsMol(:,21) = (/ 2, 14, 12,  9, 10 /)
      DihrsMol(:,22) = (/ 2, 14, 12,  9, 11 /)
      DihrsMol(:,23) = (/ 1, 15, 12,  9,  7 /)
      DihrsMol(:,24) = (/ 2, 15, 12,  9, 10 /)
      DihrsMol(:,25) = (/ 2, 15, 12,  9, 11 /)
      DihrsMol(:,26) = (/ 2, 16, 15, 12,  9 /)
      DihrsMol(:,27) = (/ 2, 16, 15, 12, 13 /)
      DihrsMol(:,28) = (/ 2, 16, 15, 12, 14 /)
      DihrsMol(:,29) = (/ 2, 17, 15, 12,  9 /)
      DihrsMol(:,30) = (/ 2, 17, 15, 12, 13 /)
      DihrsMol(:,31) = (/ 2, 17, 15, 12, 14 /)
      DihrsMol(:,32) = (/ 2, 18, 15, 12,  9 /)
      DihrsMol(:,33) = (/ 2, 18, 15, 12, 13 /)
      DihrsMol(:,34) = (/ 2, 18, 15, 12, 14 /)
      !
   case ( 2 ) ! Gamma hexene
      molform    = 'H3C-CH2-CH=CH-CH2-CH3'
      molname    = 'Hexene_gamma'
      den0       = 0.673D0
      natomMol   = 18
      nbondsMol  = 17
      nanglsMol  = 30
      ndihrsMol  = 34
      !
      allocate( bondsMol(3, nbondsMol), AnglsMol(4, nanglsMol),                &
      &         dihrsMol(5, ndihrsMol), attypMol(natomMol),                    &
      &         atposMol(3,natomMol), chargeMol(natomMol), molIDMol(natomMol) )                
      !
      attypMol = (/ 1, 4, 4, 4, 1, 4, 4, 2, 3, 2, 3, 1, 4, 4, 1, 4, 4, 4 /)
      atposMol(:, 1)= (/  -7.84520D0,   1.45030D0,  -0.00020D0 /)
      atposMol(:, 2)= (/  -8.72790D0,   2.08970D0,  -0.00060D0 /)
      atposMol(:, 3)= (/  -7.85330D0,   0.82120D0,   0.88990D0 /)
      atposMol(:, 4)= (/  -7.85310D0,   0.82050D0,  -0.88980D0 /)
      atposMol(:, 5)= (/  -6.57780D0,   2.32400D0,  -0.00030D0 /)
      atposMol(:, 6)= (/  -6.57000D0,   2.95270D0,  -0.89070D0 /)
      atposMol(:, 7)= (/  -6.57030D0,   2.95340D0,   0.88960D0 /)
      atposMol(:, 8)= (/  -5.36350D0,   1.44280D0,   0.00020D0 /)
      atposMol(:, 9)= (/  -5.47650D0,   0.35870D0,   0.00060D0 /)
      atposMol(:,10)= (/  -4.14840D0,   1.98560D0,   0.00010D0 /)
      atposMol(:,11)= (/  -4.03540D0,   3.06970D0,  -0.00030D0 /)
      atposMol(:,12)= (/  -2.93410D0,   1.10430D0,   0.00070D0 /)
      atposMol(:,13)= (/  -2.94180D0,   0.47570D0,   0.89110D0 /)
      atposMol(:,14)= (/  -2.94160D0,   0.47490D0,  -0.88920D0 /)
      atposMol(:,15)= (/  -1.66670D0,   1.97810D0,   0.00050D0 /)
      atposMol(:,16)= (/  -0.78400D0,   1.33870D0,   0.00090D0 /)
      atposMol(:,17)= (/  -1.65880D0,   2.60790D0,   0.89010D0 /)
      atposMol(:,18)= (/  -1.65860D0,   2.60720D0,  -0.88960D0 /)
      !
      chargeMol = (/ -0.180D0, 0.060D0,  0.060D0,  0.060D0, -0.120D0,  &
      &              0.060D0,  0.060D0, -0.115D0,  0.115D0, -0.115D0,  &
      &              0.115D0, -0.120D0,  0.060D0,  0.060D0, -0.180D0,  &
      &              0.060D0,  0.060D0,  0.060D0                      /) 
      !                  Type,  atom 1, atom 2
      bondsMol(:,  1) = (/  2,   2,   1/)
      bondsMol(:,  2) = (/  2,   3,   1/)
      bondsMol(:,  3) = (/  2,   4,   1/)
      bondsMol(:,  4) = (/  1,   5,   1/)
      bondsMol(:,  5) = (/  2,   6,   5/)
      bondsMol(:,  6) = (/  2,   7,   5/)
      bondsMol(:,  7) = (/  3,   8,   5/)
      bondsMol(:,  8) = (/  5,   9,   8/)
      bondsMol(:,  9) = (/  4,  10,   8/)
      bondsMol(:, 10) = (/  5,  11,  10/)
      bondsMol(:, 11) = (/  3,  12,  10/)
      bondsMol(:, 12) = (/  2,  13,  12/)
      bondsMol(:, 13) = (/  2,  14,  12/)
      bondsMol(:, 14) = (/  1,  15,  12/)
      bondsMol(:, 15) = (/  2,  16,  15/)
      bondsMol(:, 16) = (/  2,  17,  15/)
      bondsMol(:, 17) = (/  2,  18,  15/)
      !                   Type,  atom1,  atom2, atom3
      AnglsMol(:,  1) = (/  3,   2,   1,   3 /)
      AnglsMol(:,  2) = (/  3,   2,   1,   4 /)
      AnglsMol(:,  3) = (/  2,   2,   1,   5 /)
      AnglsMol(:,  4) = (/  3,   3,   1,   4 /)
      AnglsMol(:,  5) = (/  2,   3,   1,   5 /)
      AnglsMol(:,  6) = (/  2,   4,   1,   5 /)
      AnglsMol(:,  7) = (/  2,   1,   5,   6 /)
      AnglsMol(:,  8) = (/  2,   1,   5,   7 /)
      AnglsMol(:,  9) = (/  4,   1,   5,   8 /)
      AnglsMol(:, 10) = (/  3,   6,   5,   7 /)
      AnglsMol(:, 11) = (/  5,   6,   5,   8 /)
      AnglsMol(:, 12) = (/  5,   7,   5,   8 /)
      AnglsMol(:, 13) = (/  8,   5,   8,   9 /)
      AnglsMol(:, 14) = (/  6,   5,   8,  10 /)
      AnglsMol(:, 15) = (/  7,   9,   8,  10 /)
      AnglsMol(:, 16) = (/  7,   8,  10,  11 /)
      AnglsMol(:, 17) = (/  6,   8,  10,  12 /)
      AnglsMol(:, 18) = (/  8,  11,  10,  12 /)
      AnglsMol(:, 19) = (/  5,  10,  12,  13 /)
      AnglsMol(:, 20) = (/  5,  10,  12,  14 /)
      AnglsMol(:, 21) = (/  4,  10,  12,  15 /)
      AnglsMol(:, 22) = (/  3,  13,  12,  14 /)
      AnglsMol(:, 23) = (/  2,  13,  12,  15 /)
      AnglsMol(:, 24) = (/  2,  14,  12,  15 /)
      AnglsMol(:, 25) = (/  2,  12,  15,  16 /)
      AnglsMol(:, 26) = (/  2,  12,  15,  17 /)
      AnglsMol(:, 27) = (/  2,  12,  15,  18 /)
      AnglsMol(:, 28) = (/  3,  16,  15,  17 /)
      AnglsMol(:, 29) = (/  3,  16,  15,  18 /)
      AnglsMol(:, 30) = (/  3,  17,  15,  18 /)
      !                      Type, atom1, atom2, atom3, atom4
      DihrsMol(:,  1 ) = (/  2,   6,   5,   1,   2 /)
      DihrsMol(:,  2 ) = (/  2,   6,   5,   1,   3 /)
      DihrsMol(:,  3 ) = (/  2,   6,   5,   1,   4 /)
      DihrsMol(:,  4 ) = (/  2,   7,   5,   1,   2 /)
      DihrsMol(:,  5 ) = (/  2,   7,   5,   1,   3 /)
      DihrsMol(:,  6 ) = (/  2,   7,   5,   1,   4 /)
      DihrsMol(:,  7 ) = (/  3,   8,   5,   1,   2 /)
      DihrsMol(:,  8 ) = (/  3,   8,   5,   1,   3 /)
      DihrsMol(:,  9 ) = (/  3,   8,   5,   1,   4 /)
      DihrsMol(:, 10 ) = (/  4,   9,   8,   5,   1 /)
      DihrsMol(:, 11 ) = (/  4,   9,   8,   5,   6 /)
      DihrsMol(:, 12 ) = (/  4,   9,   8,   5,   7 /)
      DihrsMol(:, 13 ) = (/  5,  10,   8,   5,   1 /)
      DihrsMol(:, 14 ) = (/  6,  10,   8,   5,   6 /)
      DihrsMol(:, 15 ) = (/  6,  10,   8,   5,   7 /)
      DihrsMol(:, 16 ) = (/  7,  11,  10,   8,   5 /)
      DihrsMol(:, 17 ) = (/  7,  11,  10,   8,   9 /)
      DihrsMol(:, 18 ) = (/  7,  12,  10,   8,   5 /)
      DihrsMol(:, 19 ) = (/  7,  12,  10,   8,   9 /)
      DihrsMol(:, 20 ) = (/  6,  13,  12,  10,   8 /)
      DihrsMol(:, 21 ) = (/  4,  13,  12,  10,  11 /)
      DihrsMol(:, 22 ) = (/  6,  14,  12,  10,   8 /)
      DihrsMol(:, 23 ) = (/  4,  14,  12,  10,  11 /)
      DihrsMol(:, 24 ) = (/  5,  15,  12,  10,   8 /)
      DihrsMol(:, 25 ) = (/  4,  15,  12,  10,  11 /)
      DihrsMol(:, 26 ) = (/  3,  16,  15,  12,  10 /)
      DihrsMol(:, 27 ) = (/  2,  16,  15,  12,  13 /)
      DihrsMol(:, 28 ) = (/  2,  16,  15,  12,  14 /)
      DihrsMol(:, 29 ) = (/  3,  17,  15,  12,  10 /)
      DihrsMol(:, 30 ) = (/  2,  17,  15,  12,  13 /)
      DihrsMol(:, 31 ) = (/  2,  17,  15,  12,  14 /)
      DihrsMol(:, 32 ) = (/  3,  18,  15,  12,  10 /)
      DihrsMol(:, 33 ) = (/  2,  18,  15,  12,  13 /)
      DihrsMol(:, 34 ) = (/  2,  18,  15,  12,  14 /)
      !
   case default
      info = 'Unknown beta alkene, exit!'
      call error( subname, info, 1 )
   end select
   !
   molIDMol  = 1
   nMol      = 1
   molmass   = 0.D0
   do i = 1, natomMol
      molmass = molmass + EMass( attypMol(i) )
   enddo
   forall ( i=1:3 ) minbox(i) = maxval(atposMol(i,:)) - minval(atposMol(i,:))
   minlen = maxval( minbox )
   allocate( rotatedMol( 3, natomMol) )
   !
   subname = subnamep
return
end subroutine
